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LATTICE  BOLTZMANN  EQUATION  ON  A  2D  RECTANGULAR  GRID 

M’HAMED  BOUZIDI*,  DOMINIQUE  D’HUMIERESt,  PIERRE  LALLEMANDi,  AND  LI-SHI  LUO§ 

Abstract.  We  construct  a  multi-relaxation  lattice  Boltzmann  model  on  a  two-dimensional  rectangular 
grid.  The  model  is  partly  inspired  by  a  previous  work  of  Koelman  to  construct  a  lattice  BGK  model  on  a 
two-dimensional  rectangular  grid.  The  linearized  dispersion  equation  is  analyzed  to  obtain  the  constraints 
on  the  isotropy  of  the  transport  coefficients  and  Galilean  invariance  for  various  wave  propagations  in  the 
model.  The  linear  stability  of  the  model  is  also  studied.  The  model  is  numerically  tested  for  three  cases:  (a) 
a  vortex  moving  with  a  constant  velocity  on  a  mesh  with  periodic  boundary  conditions;  (b)  Poiseuille  flow 
with  an  arbitrary  inclined  angle  with  respect  to  the  lattice  orientation;  and  (c)  a  cylinder  asymmetrically 
placed  in  a  channel.  The  numerical  results  of  these  tests  are  compared  with  either  analytic  solutions  or  the 
results  obtained  by  other  methods.  Satisfactory  results  are  obtained  for  the  numerical  simulations. 

Key  words,  generalized  lattice  Boltzmann  equation,  rectangular  meshes,  stability  analysis,  dispersion 
equation,  Taylor  vortex,  Poiseuille  flow,  flow  past  a  cylinder  in  channel 

Subject  classification.  Fluid  Mechanics 

1.  Introduction.  Historically  originating  from  the  lattice  gas  automata  (LG A)  introduced  by  Frisch, 
Hasslacher,  and  Pomeau  [3],  the  lattice  Boltzmann  equation  (LBE)  has  recently  become  an  alternative 
method  for  computational  fluid  dynamics.  The  essential  ingredients  in  any  lattice  Boltzmann  models  which 
are  required  to  be  completely  specified  are:  (i)  a  discrete  lattice  space  on  which  fluid  particles  reside;  (ii)  a 
set  of  discrete  velocities  (often  going  from  one  node  to  its  nearest  neighbors)  to  represent  particle  advection; 
and  (iii)  a  set  of  rules  for  the  redistribution  of  particles  residing  on  a  node  to  mimic  collision  processes  in 
a  real  fluid.  Fluid-boundary  interactions  are  usually  approximated  by  simple  reflections  of  the  particles  by 
solid  interfaces. 

In  a  hydrodynamic  simulation  by  using  the  lattice  Boltzmann  equation,  one  solves  the  evolution  equations 
of  the  distribution  functions  of  fictitious  fluid  particles  colliding  and  moving  synchronously  on  a  highly 
symmetric  lattice  space.  The  highly  symmetric  lattice  space  is  a  result  of  the  discretization  of  particle 
velocity  space  and  the  condition  for  synchronous  motions.  That  is,  the  discretizations  of  time  and  particle 
phase  space  are  coherently  coupled  together.  This  makes  the  evolution  of  the  lattice  Boltzmann  equation 
very  simple  —  it  consists  in  only  two  steps:  collision  and  advection.  One  immediate  limitation  of  the 
LBE  method  is  due  to  its  use  of  highly  symmetric  regular  lattice  mesh,  which  are  usually  triangular  or 
square  lattices  in  two  dimensions  and  cubic  in  three  dimensions.  Obviously  this  is  a  serious  obstacle  to  its 
applications  in  many  areas  of  computational  fluid  dynamics.  To  deal  with  complex  computational  domains, 
various  proposals  have  been  made  to  use  grids  that  are  better  suited  to  fit  boundaries  or  to  adapt  meshes 
according  to  the  physics  of  the  system. 
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It  has  been  shown  recently  that  the  lattice  Boltzmann  equation  is  indeed  a  special  finite  difference 
form  of  the  continuous  Boltzmann  equation  with  some  drastic  approximations  tailored  for  hydrodynamic 
simulations  [4,  5,  7].  This  makes  the  lattice  Boltzmann  method  more  amenable  to  incorporate  body-fitted 
meshes  [8,  9]  or  grid  refinement  techniques  [10].  In  most  cases,  the  regular  lattice  mesh  is  abandoned  by 
decoupling  the  spatial-temporal  discretizations  and  the  discrete  velocity  set,  so  that  interpolations  can  be 
used  in  addition  to  the  advection  on  a  non-regular  or  non-uniform  mesh.  However,  interpolations  introduce 
additional  numerical  viscosities  and  other  artifacts  into  the  lattice  Boltzmann  method  [11].  Therefore,  it  is 
highly  desirable  to  construct  lattice  Boltzmann  models  with  arbitrary  mesh  and  free  of  interpolations  [2,  12]. 

In  this  paper  we  shall  consider  a  two  dimensional  model  on  a  rectangular  grid  with  an  aspect  ratio 
of  a  =  6yfSx,  where  0  <  a  <  1.  The  model  is  inspired  in  part  by  a  previous  work  of  Koelman  [2]  who 
proposed  a  general  scheme  to  construct  lattice  BGK  models  with  given  discrete  velocity  sets  based  on  a  low 
Mach  number  expansion  of  the  Maxwellian  equilibrium  distribution  function.  Conservation  and  symmetry 
constraints  are  imposed  to  fix  the  parameters  in  the  equilibrium  distribution  function.  Koelman’s  model  is 
essentially  a  variation  of  the  lattice  BGK  model  [13,  14].  As  we  shall  show,  the  transport  coefficients  of  this 
model  are  generally  anisotropic  when  a  ^  1  [15]. 

We  use  the  generalized  lattice  Boltzmann  equation  with  multiple  relaxation  times  due  to  d’Humieres  [1], 
instead  of  the  standard  lattice  BGK  model  [13,  14].  The  generalized  LBE  model  has  the  freedom  of  multiple 
relaxations  which  can  be  independent  or  coupled  together.  This  allows  one  to  optimize  the  overall  properties 
of  the  model  through  suitable  compensation  of  inadequate  behaviors.  We  shall  study  the  time  evolution 
of  plane  waves  by  analyzing  the  linearized  dispersion  equation  of  the  model  [11].  This  analysis  allows  us 
to  obtain  the  conditions  under  which  the  model  can  be  used  to  simulate  the  Navier-Stokes  equation,  i.e., 
the  model  is  Galilean  invariant  and  isotropic  up  to  a  certain  order  in  wave-number  k.  We  show  that  severe 
stability  constraints  are  due  to  Galilean  invariance  and  isotropy  of  transport  coefficients.  This  demonstrates 
the  difficulty  in  the  endeavor  of  constructing  a  lattice  Boltzmann  model  with  arbitrary  grid.  Simulations  of 
non-trivial  cases  are  presented  to  demonstrate  the  qualities  and  defects  of  the  model. 

We  organize  the  paper  as  follows.  Section  2  describes  the  proposed  model  on  a  rectangular  grid.  Section  3 
shows  a  detailed  analysis  of  the  dispersion  equation.  The  wave-number  dependence  of  Galilean  coefficient  and 
attenuation  coefficients  are  computed  explicitly  to  obtain  the  conditions  under  which  the  model  is  Galilean 
invariant  and  isotropic.  Section  4  provides  examples  of  numerical  simulations  using  the  proposed  model:  (a) 
a  vortex  moving  with  a  uniform  velocity  in  a  periodic  system;  (b)  Poiseuille  flow  with  the  boundaries  along 
arbitrary  direction  with  respect  to  the  underlying  lattice;  and  (c)  flow  past  a  cylinder  asymmetrically  placed 
in  a  channel.  Section  5  concludes  the  paper. 


2.  Definition  of  the  model.  We  consider  a  two-dimensional  LBE  model  with  nine  discrete  velocities 
(the  D2Q9  model)  on  a  rectangular  grid  of  dimensions  1  and  a.  (In  what  follows  all  quantities  are  given  in 
non-dimensional  units,  normalized  by  the  lattice  unit  6x-)  In  the  advection  step  of  the  lattice  Boltzmann 
equation,  particles  move  from  one  node  of  the  grid  to  one  of  its  neighbors  as  illustrated  in  Eig.  1.  The 
discrete  velocities  are  given  by 


ea=  < 


(0,  0) , 

(cos[(q:  —  l)7r/2],  a  sin[(Q:  —  l)7r/2]), 
(cos[(2q:  —  9)7r/4],  a  sin[(2Q:  —  9)7r/4])'\/2, 


O'  =  0, 

a  =  1-4, 
a  =  5-8, 


(2.1) 
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where  the  duration  of  the  time  step  6t  is  assumed  to  be  unity.  At  any  time  the  LBE  fluid  is  then 
characterized  by  the  populations  of  the  nine  velocities  at  each  node  of  the  computational  domain 

|/(rj,  tn))  =  {Mvj,  t„),  tn),  ■■■  ,  faivj,  t„)y ,  (2.2) 

where  T  is  the  transpose  operator.  Here  upon  the  Dirac  notation  of  bra,  (-j,  and  ket,  j-),  vectors  is  used  to 
denote  row  column  and  row  vectors,  respectively.  The  time  evolution  of  the  state  of  the  fluid  follows  the 
general  equation 

\f{rj  +  e„,  tn  +  1))  =  \f{rj,  tn))  +  Mfirj,  tn))) ,  (2.3) 


where  collisions  are  symbolically  represented  by  the  operator  Q. 


r - ^ - i 
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Fig.  1.  Discrete  velocities  of  the  nine-velocity  model  on  a  rectangular  grid.  The  aspect  ratio  of  the  rectangle  is  SyjSx  =  a. 


We  shall  use  the  generalized  lattice  Boltzmann  equation  introduced  by  d’Humieres  [1],  in  which  the 
collision  process  is  executed  in  moment  space  M.  The  mapping  between  moment  space  and  discrete  velocity 
space  V  spanned  by  {eQ,}  is  one-to-one  and  defined  by  the  linear  transformation  M  which  maps  a  vector  |/) 
in  V  to  a  vector  |/)  in  M,  i.e., 

|/)  =  M|/),  and  |/)  =  M-i|/).  (2.4) 


To  reflect  the  underlying  symmetries  appearing  in  both  the  Chapman-Enskog  expansion  and  the  dispersion 
equation,  M  is  constructed  as  the  following 
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(2.5) 


where  -I- 1,  <^2  =  1  —  2a^,  ips  =  —  2,  ip4  =  —  1,  <^5=0^-!-  2,  and  ipe  =  —(1  -I-  2a^). 

The  components  of  the  row  vector  {mp\  in  matrix  M  are  polynomials  of  the  x  and  y  components  of  the 
velocities  {eQ,},  ea,x  and  6^,^.  The  vectors  {mp\,  /3  =  1,  2,  •  •  •  ,  9,  are  orthogonalized  by  the  Gram-Schmidt 
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procedure  in  a  carefully  considered  order.  The  first  three  orthogonal  vectors  correspond  to  the  mass,  x-  and 
j/-momentum  modes:  (mi|  =  (lie-in,  {mi\  =  {ea,x\^  and  (mel  =  {ea,y\-  The  above  expressions  prescribe  the 
components  of  (mi  | ,  (m4 1 ,  and  (me  | .  These  three  vectors  span  the  hydrodynamic  subspace  of  the  eigen-space 
of  the  collision  operator  for  a  two-dimensional  athermal  LBE  model.  The  remaining  six  vectors  span  the 
kinetic  subspace.  The  vector  (m2|  =  (3||ea|p  —  2(1  -|-  a^)||ea||°|  is  constructed  by  orthogonalizing  the  energy 
mode  (||ea|p|.  Similarly,  vectors  (msl  and  (myl  are  respectively  built  upon  (eQ,_3;||eQ,|p|  and  (eQ,,j,||eQ,|p|; 
(mgl  is  constructed  upon  (e^  y|  and  (mgj  =  {ea,xSa,y\i  and  finally  {mz\  is  obtained  from  (||ea||'^|.  By 

means  of  their  construction,  the  row  vectors  in  M  are  mutually  orthogonal,  but  they  are  not  normalized, 
their  norms  being  chosen  to  simplify  algebraic  manipulations.  When  a  =  1,  M  reduces  to  that  for  the  D2Q9 
model  on  a  square  grid  with  a  different  normalization  of  \pxx)  [H]-  Therefore  the  proposed  model  can  be 
considered  as  a  generalization  of  the  model  on  a  square  lattice.  It  should  be  noted  that  when  a  1,  there 
are  three  nonzero  (kinetic)  energy  levels  in  the  model  which  introduce  additional  degrees  of  freedom  into 
the  model  and  extra  care  must  be  taken  in  the  construction  of  (mg],  (mg],  and  (rngj,  i.e.,  they  must  be 
orthogonalized  with  the  Gram-Schmidt  procedure  in  the  particular  order  of  (m2 1 ,  (mg  | ,  and  (m3 1 . 

It  is  interesting  to  note  that  the  moments  {mp\f)  have  a  physical  interpretation.  The  matrix  M  so 
constructed  in  the  above  naturally  leads  to  the  moment  vector  in  moment  space  M  as  the  following: 


I/)  —  (Pj  6,£,  jx,  Qxi  jy7  Qyi  Pxx7  Pxy)  j 


(2.6) 


where  p  is  the  density,  e  is  related  to  the  kinetic  energy,  e  is  related  to  the  kinetic  energy  squared  for  a  =  1 
(but  has  no  obvious  physical  meaning  when  a  1),  jx  and  jy  are  x  and  y  components  of  the  momentum 
density,  Qx  and  qy  are  proportional  to  the  x  and  y  components  of  the  energy  flux,  and  Pxx  and  Pxy  are 
proportional  to  the  diagonal  and  off-diagonal  components  of  the  viscous  stress  tensor. 

For  the  collision  process,  we  propose  to  use  the  following  equilibrium  distribution  functions  of  the 
(nonconserved)  moments,  which  depend  only  on  the  conserved  moments,  i.e.,  p,  and  j  =  {jx,  jy). 


=,(eq) 
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i(eq) 
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2{3cl  -  1  -  a^)p+^{jl  +jy) , 
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2^2jy  > 

[3(„=  +  l)cl-  2„=]  ,  +  2  („=,=  - 

1  .  . 

'^JxJy  j 


(2.7a) 

(2.7b) 

(2.7c) 

(2.7d) 

(2.7e) 

(2.7f) 


where  the  coupling  coefficient  between  p^^^^  and  p  (which  vanishes  in  the  standard  D2Q9  LBE  model)  is 
introduced  to  obtain  the  isotropy  of  the  sound  speed.  The  values  of  the  coupling  constants  (0:3,  ci,  and  C2) 
in  the  above  equilibria  are  obtained  by  optimizing  isotropy  and  and  stability  of  the  model  [11].  It  should  be 
noted  that  the  energy  is  not  considered  as  a  conserved  quantity  here  because  the  model  is  athermal.  (The 
model  does  not  possess  sufficient  degrees  of  freedom  to  accommodate  the  dynamics  of  locally  isotropic  heat 
transport). 

In  what  follows  the  idea  of  the  “incompressible”  LBE  [6]  is  applied  to  the  above  equilibria  so  that  p 
is  replaced  by  a  constant  po  in  the  denominators  of  equations  (2.7a),  (2.7e),  and  (2.7f).  This  choice  allows 
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for  better  comparison  with  other  incompressible  simulations  and  simpler  algebra  while  retaining  correct 
acoustics. 

The  collision  process  is  modeled  by  the  following  relaxation  equations 

ir)  =  l/)-S[|/)-|/^'^"')],  (2.8) 

where  |/*)  denotes  the  post-collision  state,  and  S  is  the  diagonal  relaxation  matrix 

S  =  diag(0,  S2,  S3,  0,  S5,  0,  S7,  ss,  sg)  ■  (2.9) 

The  model  reduces  to  the  usual  lattice  BGK  model  if  all  the  relaxation  parameters  are  set  to  be  a  single 
relaxation  time  r  (and  a  =  1),  i.e.,  Sa  =  1/r.  It  should  be  stressed  that  the  relaxation  parameters  are  not 
independent,  as  shown  in  the  next  section.  The  constraints  of  isotropy  lead  to  the  coupling  between  these 
relaxation  parameters  [11].  Obviously,  the  usual  lattice  BGK  model  does  not  possess  the  freedom  for  such 
couplings,  therefore  it  would  not  work  on  a  rectangular  grid. 

3.  Analysis  of  linearized  dispersion  eqnation.  The  analysis  presented  in  what  follows  is  similar 
to  that  presented  in  Ref.  [11]  where  the  goal  of  the  work  was  to  determine  the  stability  conditions  for  the 
coupling  coefficients  Q!2  and  as,  and  the  constraints  on  the  relaxation  parameters  Sa- 

We  consider  a  system  of  size  iVj,  x  Ny  with  periodic  boundary  conditions  and  look  for  small  amplitude 
solutions  in  the  presence  of  a  uniform  flow  [for  given  values  of  p  and  V  =  (I4,  Vy)  =  J/p].  For  a  wave  vector 
k  in  the  reciprocal  space  of  the  computational  domain,  we  search  for  solutions 

fa{r,t)  on  exp{—ik  ■  r  +  zt) .  (3.1) 

To  first  order  in  fc,  we  have  the  following  linearized  dispersion  equation: 

det(K(i)+M-iCM-.2l)  =  0, 

where  I  is  the  identity  operator,  is  the  linearized  advection  operator  which 

=  diag(0,  ik  ■  ei,  . . . ,  ik  ■  eg) , 

and  C  is  the  linearized  collision  operator 
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where 


(3.2) 

is  a  diagonal  matrix 

(3.3) 


ai  =  2[Zcl  -  (1  +  a^)]  , 

h^[2a^-3c^(a^  +  l)]  . 


Oih  = 


5 


(3.5a) 

(3.5b) 


The  linearized  dispersion  equation  (3.2)  can  be  solved  by  perturbation  technique  in  power  series  of  k  [11]. 
To  ensure  isotropy  and  Galilean  invariance  in  the  limit  of  fc  ^  0,  we  need  to  solve  the  linearized  dispersion 
equation  up  to  . 

In  the  first  order  of  fe,  three  solutions  are  obtained:  one  corresponds  to  transverse  excitations  which  are 
convected  with  the  uniform  speed  of  the  fluid  k-V  jk,  whereas  the  other  two  are  acoustic  waves  with  phase 
velocity  ic^,  where  the  speed  of  sound  Cg  can  be  chosen  within  limits  that  is  deferred  to  later  discussion. 
The  sound  waves  also  have  the  correct  dependence  on  the  applied  uniform  velocity  V  of  the  fluid  up  to  the 
first  order  in  V,  i.e.,  Cg  ^  Cg  ±V  cos^,  where  ^  is  the  angle  between  k  and  V.  The  nonlinear  terms  in  the 
equilibria  of  Eqs.  (2.7a)  -  (2.7f)  provide  the  correct  Galilean  coefficients  for  both  transverse  and  longitudinal 
waves. 

In  the  second  order  (in  k)  of  the  solutions  of  the  dispersion  equation,  the  constraints  on  the  isotropy  of 
the  transport  coefficients  for  the  hydrodynamic  modes  lead  to 


C2  +  4(1  -  a^) 

Cl  — - 5 - 

and  the  following  relationships  between  the  relaxation  parameters 

j_  _ _ 2(4+C2)[(12c^-C2)(l  +  a^)-2(5a=^  +  2)] _ j_ 

§2  (l  +  a2)(l  +  C2-3a2)(c2  +  10-12c2)  +  6[a"‘(c2-2)-3(a2-l)]  gg  ’ 

1  _  2(4  +  C2)[(12c^ -C2)(l  +  a^)-2(3a"‘  +  5a^  +  5)]  1 

sg  (l  +  a2)(l  +  c2-3a2)(c2  +  10-12c2)  +  6[a'‘(c2-2)-3(a2-l)]  sq  ’ 


(3.6) 


(3.7a) 

(3.7b) 


where  1/sa  =  (1/sa  —  1/2).  The  coupling  between  S2  and  sg  is  required  only  when  a  ^  1.  The  kinematic 
shear  viscosity  v  and  the  kinematic  bulk  viscosity  (  are 


4  +  C2 


C  =  ^(7  +  3a^  +  C2-12cl)  - 


12 


S2 


(3.8a) 

(3.8b) 


For  a  given  a,  the  speed  of  sound  and  C2  must  be  chosen  such  that  the  shear  and  bulk  viscosities  are  positive 
and  the  Eqs.  (3.7a)  and  (3.7b)  lead  to  positive  values  for  S2  and  sg. 

The  values  of  Cg  and  C2  which  optimize  the  isotropy  and  stability  of  the  model,  depending  on  the  grid 
aspect  ratio  a,  are  determined  by  the  linear  analysis  of  the  model  [11].  In  the  case  of  square  grid,  i.e., 
a  =  1,  we  have  found  =  1/3  and  C2  =  —2.  This  result  coincides  with  the  one  given  in  Ref.  [11]  and  the 
relationship  between  sg  and  sg  given  by  Eq.  (3.7b),  and  the  shear  and  bulk  viscosities  given  by  Eqs.  (3.8a) 
and  (3.8b),  all  reduces  to  the  previous  results  for  a  square  lattice  where  ci  =  —2  [see  Eqs.  (40),  (41),  (42), 
and  (43)  in  Ref.  [11]  for  Cg,  sg(sg),  i/,  and  respectively].  However,  the  coupling  between  S2  and  sg  is  unique 
to  the  model  on  a  rectangular  grid.  This  coupling  is  due  to  the  dependence  of  on  p,  which  in  turn  leads 
the  term  0:5 sg  in  the  linearized  collision  operator  C  in  Eq.  (3.4).  Finally  note  that  0:3  has  little  influence  and 
is  set  to  be  equal  to  —2. 

The  linearized  dispersion  equation  can  be  solved  numerically  for  any  value  of  k  to  determine  the  linear 
stability  of  the  system  by  computing  the  rate  of  growth  of  spatially  periodic  excitations  superimposed  to 
a  uniform  flow  of  constant  velocity  V,  as  previously  shown  in  the  case  of  a  square  grid  [11].  Through  this 
analysis  it  is  found  that  the  present  model  is  much  less  stable  than  the  square  one,  i.e.,  the  stable  region  in 
parameter  space  of  V  and  is  much  smaller  than  that  for  the  model  with  a  square  lattice.  For  instance, 
when  a  =  1/2,  a  stability  condition  is  that  V  <  0.05,  whereas  for  the  model  with  a  square  lattice  (a  =  1),  the 
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same  stability  condition  is  that  V  <  0.20.  One  reason  for  this  is  that  in  general  the  sound  speed  Cg  decreases 
with  the  aspect  ratio  a;  for  instance,  when  a  =  1/2  the  optimal  speed  of  sound  is  about  0.377,  which  is 
different  from  the  usual  Cg  =  1/V3  ~  0.577  on  a  square  lattice.  Therefore,  the  local  velocity  magnitude 
must  be  decreased  accordingly  to  keep  the  local  Mach  in  check  so  that  the  low  Mach  number  approximation 
remains  valid.  This  means  that  the  present  model  will  have  limited  ability  to  simulate  flows  even  at  moderate 
Reynolds  numbers.  In  addition  when  using  a  combination  of  rectangular  and  square  grids  (the  simplest  case 
of  grid  refinement  in  one  direction)  in  the  situation  where  acoustic  propagation  is  important,  it  will  not  be 
possible  to  choose  an  optimal  value  of  the  sound  speed  for  the  two  different  grids. 

We  would  like  to  note  that,  although  there  is  no  simple  interpretation  of  the  instability  of  the  LBE 
models  due  to  the  presence  of  a  uniform  velocity  V,  information  on  the  instability  can  be  obtained  by 
analyzing  the  velocity  dependence  of  the  attenuation  of  sound  waves  using  the  linearized  dispersion  equation 
[11]. 

Let  us  consider  the  case  where  the  uniform  velocity  is  parallel  to  the  wave  vector  k  with  a  polar  angle 
(between  k  and  a;-axis).  For  small  values  of  k  and  the  particular  choice  of  C2 

C2  =  (a2-3),  (3.9) 


we  have  the  following  results.  The  transverse  mode  has  phase  velocity  v±  =  V  and  its  attenuation  is  given 

by 

9(1  —  a^)^  sin^  29 


7-L  =  ^ 


1  +  _y2 


sg 


6 


2[l  -  12,0?  +  +  Q(l  +  a?)cl] 


(3.10) 


For  the  longitudinal  modes,  we  obtain  as  phase  velocity  uy  =  and  attenuation  coefficient  yy  = 

(jb  +  7-l)/2,  with  (to  the  first  order  in  V) 


1  +  —  3c„ 


V 


3  Cs[l  +  8a^  +  —  12(1  +  a^)c^] 

-3  (7  +  120^  +  3a^)  +12(1  -  a^)cl  cos^9  [2  +  (1  -  a^)(2  -  3  cos^  ^')] }  j  ■ 


{(l  +  a2)(7a2  +36cj) 


(3.11) 


Contrary  to  the  case  of  square  grid,  it  is  not  possible  in  general  with  a  given  value  of  a  1  to  find  a  value  of 
Cg  for  which  the  linear  dependence  of  the  attenuation  of  acoustic  waves  on  V  can  be  eliminated  (for  a  =  1, 
this  can  be  accomplished  by  setting  =  1/3).  This  is  a  possible  cause  of  instability  in  the  model. 

4.  Simulations.  We  use  the  two-dimensional  multi-relaxation  LBF  model  on  either  a  square  grid  or 
a  rectangular  grid  for  the  following  simulations.  The  central  routine  (collision  and  advection)  is  quite  close 
to  that  for  the  standard  square  LBE  and  leads  to  similar  performances  (using  a  workstation  with  a  500 
MHz  EV6  processor,  the  overall  computation  time  per  node  and  per  time  step  is  in  the  range  0.2  to  0.4 
microsecond  depending  whether  the  cache  is  large  enough  or  not). 

4.1.  a  vortex  traveling  with  a  constant  velocity.  To  test  the  ability  of  the  present  LBE  scheme  to 
simulate  a  viscous  flow,  we  consider  the  particular  case  of  a  simple  vortex  superimposed  to  a  uniform  flow 
of  velocity  V.  We  take  as  initial  condition  for  the  flow 


uo{r,  t  =  0)  =  V  +  {yo  -  y,  X  -  xo)  too  exp[-(r  -  ro)‘^/R^] . 


(4.1) 


where  ro  =  {xq,  yo)  is  the  initial  position  of  the  vortex  center,  and  ojq  and  R  characterize  respectively  the 
amplitude  and  the  extent  of  the  vortex.  The  evolution  of  the  corresponding  macroscopic  flow  is  fairly  simple: 
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the  center  of  the  vortex  travels  with  the  velocity  V  and  the  maximum  value  of  the  vorticity  (at  the  center 
ro  +  Vt)  decays  in  time  as 

UJvnaA  ^  -  (7^2  +  ~  (I  +  4t*)2  ’  ^ 

where  t*  =  vtlR"^  is  the  dimensionless  time. 


Fig.  2.  LBE  simulation  of  a  moving  vortex.  Decay  of  the  vorticity  maximum.  The  grid  aspect  ration  a  =  lj2.  Symbol  + 
and  X  are  simulation  results  for  y  =  0  and  V  =  0.05,  respectively.  The  solid  lines  are  fitting  of  the  data  according  to  Eq.  (4.2) 
with  the  viscosity  of  value  v  =  0.9876i^o  and  v  =  0.8966i^o;  respectively,  where  vq  is  given  by  Eq.  (3.8a). 

The  system  size  is  x  Ny  =  109  x  109,  with  a  grid  aspect  ratio  a  =  1/2.  The  size  of  the  vortex  is 
R  =  6.  Values  of  other  parameters  are:  0:2  =  —3.5,  0:3  =  2.0,  C2  =  —2.9,  and  sg  =  1.8,  i.e.,  v  =  0.01018 
according  to  Eq.  (3.8a).  The  results  obtained  by  the  LBE  simulations  with  various  conditions  agree  very  well 
with  the  analytic  solution  of  the  flow  for  V  =  0.  However  when  V  increases  there  are  departures  from  the 
simple  result  of  Eq.  (4.2)  due  to  the  dependence  of  the  transport  coefficients  and  g-factor  on  V,  as  discussed 
for  the  square  grid  in  Ref.  [11].  An  example  of  such  behavior  is  demonstrated  in  Eig.  2.  Eigure  2  shows 
two  LBE  simulation  results  of  Wmax  as  a  function  of  dimensionless  time  t*  =  vtlR"^,  with  V  =  (0,  0)  and 

V  =  (0.05,  0).  Equation  (4.2)  is  used  to  fit  the  data  to  obtain  the  viscosity.  The  results  are  v  =  0.9876i/o 
and  V  =  0.8966i/o  for  I4  =  0  and  I4  =  0.05,  respectively,  where  i/q  is  given  by  Eq.  (3.8a).  There  are  two 
factors  that  attribute  to  the  correction  in  the  viscosity:  the  wave-number  fe-dependence  and  V-dependence 
of  the  transport  coefficients  [11].  The  same  simulations  are  performed  on  a  square  grid  and  the  results  are: 

V  =  0.9866i/o  and  v  =  0.8745j/o  for  =  0  and  14  =  0.05,  respectively.  It  should  be  noted  that  in  the  LBE 
simulations,  initial  conditions  include  not  only  the  conserved  quantities  such  as  the  density  and  velocity 
fields,  but  also  all  the  nonconserved  quantities  such  as  fluxes  and  the  stress,  which  can  be  obtained  from  the 
initial  velocity  field  through  a  Chapman- Enskog  analysis  of  the  model. 

4.2.  Poiseuille  flow  with  arbitrarily  inclined  walls.  The  second  test  is  the  two-dimensional  Poiseuille 
flow  with  arbitrarily  inclined  walls.  This  situation  allows  us  to  test  the  no-slip  boundary  conditions  in  the 
LBE  model.  We  consider  a  system  of  size  x  Ny  with  periodic  boundary  conditions.  The  boundaries  of 
the  channel  are  placed  with  an  arbitrary  inclined  angle  6  with  respect  to  a;-axis,  as  illustrated  in  Eig.  3. 


The  no-slip  boundary  conditions  used  here  for  the  channel  walls  are  the  interpolated  bounce-back  boundary 
conditions  proposed  in  Ref.  [16].  The  interpolated  bounce-back  boundary  conditions  combine  interpolation 
and  bounce-back  schemes  to  deal  with  boundaries  which  are  off  the  lattice  points. 


Fig.  3.  2D  Poiseuille  flow  with  arbitrary  inclined  walls.  The  system  size  is  assumed  to  be  x  Ny.  The  dises  are  grid 
points.  The  solid  lines  are  the  adveetion  lines  of  the  diserete  veloeities.  The  dashed  lines  are  the  boundaries  of  the  channel. 
The  width  of  the  ehannel  is  Ny .  The  no-slip  boundary  conditions  are  enforced  at  the  intersections  the  dashed  lines  and  the 
thin  solid  lines. 


We  first  studied  the  time  evolution  of  the  flow  starting  at  rest,  and  compared  the  results  obtained  by 
using  the  rectangular  and  square  grids.  The  time  evolution  of  velocity  fields  of  the  two  systems  agree  very  well 
with  each  other.  We  also  studied  the  momentum  transfer  at  the  boundary.  We  found  an  excellent  agreement 
between  its  measurements  for  the  square  and  the  rectangular  grids,  and  its  expected  value:  pj/I/3_LVj|,  where 
L  is  the  length  of  the  boundary,  and  9_LVj|  is  the  normal  derivative  of  the  shear  velocity  with  respect  to  the 
wall,  computed  at  the  wall. 

Note  that  when  we  compute  the  momentum  transfer  for  the  rectangular  grid,  the  components  of  the 
usual  momentum  transfer  have  to  be  multiplied  by  a  factor  a  to  account  for  the  surface  of  the  elementary 
cell  (assuming  that  all  results  are  in  non-dimensional  units  defined  on  the  square  grid).  In  order  to  better 
understand  the  origin  of  this  factor,  one  has  to  remember  that  p  and  j  are  the  mass  and  momentum  densities 
(mass  and  momentum  per  unit  surface),  while  the  momentum  transfer  has  to  be  computed  from  momentum: 
(momentum  density)  x  (cell  surface).  Usually  a  unique  regular  grid  is  used  and  the  cell  volume  can  be  taken 
as  unit  volume.  Here  however  the  surface  of  the  cells  is  equal  to  the  chosen  aspect  ratio  a  once  the  square 
cell  has  been  taken  as  unit  surface.  Indeed  this  remark  applies  to  the  next  section  when  computing  the  drag 
and  lift  coefficients. 

4.3.  flow  past  a  cylinder  asymmetrically  placed  in  a  channel.  The  third  test  we  did  was  a  two- 
dimensional  flow  past  a  cylinder  asymmetrically  placed  in  a  channel.  This  flow  has  been  used  as  a  standard 
benchmark  test  in  CFD  [17].  The  flow  configuration  is  as  follows:  a  cylinder  of  diameter  d  is  placed  in  a 
channel  of  width  4.  Id  and  length  22d,  the  center  of  the  cylinder  is  asymmetrically  (with  respect  to  the  center 
line  of  the  channel)  located  at  horizontally  2d  from  the  entrance,  and  vertically  2d  from  the  lower  wall  of  the 
channel,  as  shown  in  Fig.  4.  The  boundary  condition  at  the  entrance  is  a  Poiseuille  profile  with  average  speed 
U.  The  boundary  condition  at  the  exit  is  free  exit  with  a  total  flux  equal  to  the  input  flux.  The  bounce-back 
boundary  conditions  are  used  for  the  channel  walls,  and  the  interpolated  bounce-back  boundary  conditions 
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Table  4.1 

2D  flow  past  a  cylinder  asymmetrically  placed  in  a  channel  at  Re  =  100. 


a 

Cs 

C2 

Nx  X  Ny 

St 

/^max 

^max 

^min 

AP 

1.00 

1/V3 

-2 

709  X  132 

0.3021 

3.153 

0.926 

-1.018 

2.50 

0.85 

0.6141 

-0.80 

709  X  155 

0.3018 

3.186 

0.984 

-1.051 

2.51 

0.80 

0.5829 

-0.90 

709  X  165 

0.3020 

3.174 

0.950 

-1.062 

2.51 

0.75 

0.5412 

-1.10 

709  X  176 

0.3018 

3.173 

0.965 

-1.053 

2.51 

0.70 

0.5113 

-1.50 

709  X  188 

0.3007 

3.195 

1.013 

-1.071 

2.51 

0.65 

0.4761 

-1.70 

709  X  203 

0.3009 

3.184 

0.999 

-1.062 

2.47 

0.60 

0.4417 

-2.00 

709  X  220 

0.3009 

3.176 

1.002 

-1.053 

2.45 

0.55 

0.4086 

-2.25 

709  X  240 

0.3015 

3.189 

1.005 

-1.052 

2.42 

0.50 

0.3770 

-2.55 

709  X  264 

0.3007 

3.199 

1.019 

-1.084 

2.45 

0.45 

0.2977 

-2.90 

709  X  293 

0.2992 

3.204 

1.053 

-1.107 

2.50 

CFD  lower  bound  in  Ref.  [17] 

0.2950 

3.22 

0.99 

— 

2.46 

CFD 

upper  bound  in 

Ref.  [17] 

0.3050 

3.24 

1.01 

— 

2.50 

with  a  second  order  interpolation  [16]  are  used  for  the  boundary  of  the  cylinder.  The  Reynolds  number  for 
the  flow  is 


We  use  the  LBE  model  to  simulate  the  flow  at  Re  =  100  for  which  there  is  periodic  vortex  shedding  behind 
the  cylinder. 


22d 

(ti,  ti)  =  (0,  0) 


flow  direction  ► 


(u,  ti)  =  (0,  0) 


Fig.  4.  Configuration  of  a  2D  flow  past  a  cylinder  asymmetrically  placed  in  a  channel. 

The  flow  was  computed  on  rectangular  grids  with  several  different  values  of  the  grid  aspect  ratio  a,  and 
compared  to  the  results  with  a  square  grid.  The  measured  quantities  are  Strouhal  number  St,  maximum  drag 
(omax^  maximum  lift  coefficient  minimum  lift  coefficient  C™'",  and  the  pressure  difference  AP.  The 

results  are  summarized  in  Table  4.1.  Table  4.1  also  shows  the  lower  and  upper  bounds  of  St,  C™’',  and 

AP,  obtained  by  a  number  of  conventional  CFD  methods  presented  in  Ref.  [17].  Overall  the  LBE  simulation 
results  with  square  or  rectangular  grids  agree  well  with  each  other,  and  with  the  CED  results  in  Ref.  [17]. 
Figures  5  show  the  contours  of  the  stream  function  tp{x,  y)  and  the  vorticity  uj{x,  y)  of  the  simulations  on 
a  square  grid  of  size  x  Ny  =  1401  x  308  and  on  a  rectangular  grid  of  size  x  Ny  =  1401  x  616.  The 
relative  P^-norm  difference  of  the  two  velocity  fields  is  about  2.2  x  10“'^.  Note  that  the  aspect  ratio  for  this 
particular  calculation  is  slightly  different  from  that  shown  in  figure  4,  but  this  has  negligible  effect  for  the 
present  purpose  of  comparing  results  on  the  square  and  the  rectangular  grids. 

The  relative  fluctuation  of  Strouhal  number  St  is  well  under  1%  and  the  values  of  St  are  well  within  the 
bounds  in  Ref.  [17].  The  fluctuation  of  is  also  under  1%  but  the  values  of  are  all  slightly  lower 
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than  the  results  in  Ref.  [17].  The  fluctuation  of  AP  is  about  1%  and  the  values  of  AP  agree  well  with  the 
results  in  Ref.  [17].  The  values  of  lift  coefficient  obtained  by  the  LBE  simulations  have  a  variation  about 
±6%,  which  is  much  greater  than  the  variations  in  other  measured  quantities. 


0  200  400  600  800  1000  1200 

Fig.  5.  2D  flow  past  a  cylinder  asymmetrically  placed  in  a  channel  at  Re=  100.  Top  and  bottom  figure  show  contours  of 
the  stream  function  ip{x,  y)  and  the  vorticity  oj{x,  y)  of  the  flow,  respectively.  The  dashed  lines  are  the  simulation  results  on 
a  square  grid  of  size  Nx  x  Ny  =  1401  x  308,  and  the  solid  lines  are  that  on  a  rectangular  grid  of  size  Nx  x  Ny  =  1401  x  616. 

A  possible  origin  of  the  discrepancy  in  the  lift  coefficients  is  the  following.  The  LBE  method  is  intrin¬ 
sically  a  compressible  scheme  and  acoustic  waves  may  be  generated  by,  e.g.,  initial  conditions  that  do  not 
include  a  proper  pressure  field  or  the  flow  itself  that  generates  an  oscillating  pressure  field  as  is  the  case 
considered  here.  Eor  a  given  value  of  the  sound  speed  and  a  given  choice  of  the  boundary  conditions  at  the 
entrance  and  exit  of  the  channel  the  frequency  of  some  of  the  longitudinal  acoustic  modes  can  be  close  to 
multiples  of  the  Strouhal  frequency  in  the  flow.  This  causes  resonances  between  some  of  the  acoustic  waves 
and  the  periodic  shedding  of  vortices  by  the  cylinder.  The  coupling  between  acoustic  waves  and  vortex  shed¬ 
ding  indeed  affects  the  hydrodynamic  fields,  and  in  turn,  various  measured  quantities.  Among  the  measured 
quantities,  the  lift  coefficients  are  most  sensitive  to  this  effect.  The  mean  drag  coefficient  is  also  affected  but 
to  a  much  smaller  extent.  This  problem  is  of  broad  interest.  However  it  will  be  easier  to  study  it  with  the 
model  of  square  grid  for  which  the  speed  of  sound  and  the  bulk  viscosity  can  be  chosen  in  a  broader  range 
than  for  the  model  of  rectangular  grid.  A  detailed  study  is  beyond  the  scope  of  the  present  work  and  will 
be  addressed  elsewhere. 

5.  Conclusion  and  discussion.  In  this  paper  we  have  successfully  proposed  a  two-dimensional  nine- 
velocity  generalized  lattice  Boltzmann  model  with  multiple  relaxations  on  a  rectangular  grid  with  arbitrary 
aspect  ratio  a  =  SyjSx-  We  have  numerically  validated  the  model  by  using  the  model  to  simulate  several 
benchmark  problems,  and  have  obtained  satisfactory  results.  In  contrast  to  the  previous  two-dimensional, 
nine-velocity,  multi-relaxation  model  on  a  square  grid  [11],  the  model  on  a  rectangular  grid  is  more  prone 
to  instability,  and  the  admissible  maximum  value  of  local  velocity  magnitude  is  much  less  than  that  in 
the  model  on  a  square  grid.  It  should  also  be  stressed  that,  although  this  work  is  in  part  motivated  by  a 
previous  work  [2],  it  is  realized  that  the  nine- velocity  lattice  BGK  equation  cannot  possibly  work  properly 
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on  a  rectangular  grid.  Specifically,  the  lattice  BGK  equation  does  not  have  sufficient  degrees  of  freedom 
to  satisfy  the  constraints  imposed  by  isotropy  and  Galilean  invariance.  With  nine  discrete  velocities  in 
two-dimensions,  it  is  necessary  to  use  the  multi-relaxations  to  construct  an  LBE  model  on  a  rectangular 
grid. 

This  work  is  our  first  attempt  to  construct  a  lattice  Boltzmann  model  on  an  arbitrary  unstructured  grid. 
As  discussed  in  Ref.  [12],  one  difficulty  encountered  in  the  LBE  model  on  an  unstructured  grid  is  due  to 
the  fact  that  Vca/  e^V/  because  the  discrete  velocity  set  {e^}  has  spatial  dependence.  In  this  work,  we 
found  that  there  are  additional  issues  in  the  LBE  model  on  an  unstructured  grid  needed  to  be  addressed. 

Eirst,  we  found  that  the  local  grid  structure  severely  affects  the  local  sound  speed.  If  the  sound  speed 
varies  spatially  depending  on  local  grid  structure,  then  the  model  is  unphysical.  Correct  acoustic  propagation 
is  an  essential  part  of  the  lattice  Boltzmann  method.  Secondly,  the  constraints  of  isotropy  and  Galilean 
invariance  are  difficult  to  satisfy  by  using  the  lattice  BGK  model,  as  proposed  in  Ref.  [12],  unless  the  discrete 
velocity  set  includes  a  large  number  of  velocities.  Thirdly,  the  numerical  stability  is  severely  affected  by  the 
local  grid  structure  even  for  uniform  structured  grid,  as  we  have  demonstrated  in  this  work.  Stability  is  of 
key  importance  to  an  effective  lattice  Boltzmann  algorithm.  However,  we  have  not  yet  developed  a  method 
to  systematically  improve  the  stability  of  the  lattice  Boltzmann  method.  We  believe  that  the  aforementioned 
issues  must  be  resolved  before  we  can  construct  a  lattice  Boltzmann  model  on  an  arbitrary  unstructured 
grid. 
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